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We report on the direct numerical measurements of the conductivity of graphene monolayer. Our 
numerical simulations are performed in the effective lattice field theory with noncompact 3+1- 
dimensional Abelian lattice gauge fields and 2 + 1-dimensional staggered lattice fermions. The 
conductivity is obtained from the Green-Kubo relations using the Maximum Entropy Method. We 
find that in a phase with spontaneously broken sublattice symmetry the conductivity rapidly de- 
creases. For the largest value of the coupling constant used in our simulations g — 4.5, the DC 
conductivity is less than the DC conductivity in the weak-coupling phase (at g < 3.5) by at least 
three orders of magnitude. 
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1. INTRODUCTION 

Graphene, a single layer of carbon atoms which form a 
two-dimensional honeycomb lattice, is probably the most 
widely discussed material in modern condensed matter 
physics. A peculiar feature of charge carriers in graphene 
is that their energy spectrum near the Fermi point is 
similar to that of the free 2+1-dimcnsional masslcss Dirac 
fermions. This explains the unusual transport properties 
of graphene, such as Klein tunneling or novel types of the 
quantum Hall effect [1-3]. 

The four spinor components of these Dirac fermions 
correspond to charge carriers which are localized on one 
of the two elementary rhombic sublatticcs of the honey- 
comb lattice and which are close to one of the two dis- 
tinct Fermi points in the Brillouin zone of graphene. In 
this low-energy description two components of the non- 
relativistic spin are treated as two independent fermionic 
flavors. Interactions between fermions are mediated by 
electromagnetic fields which propagate freely in 3 + 1- 
dimcnsional space. The strength of electromagnetic in- 
teractions can be controlled by placing graphene layers 
on substrates with different dielectric permittivities. 

Since charge carriers in graphene propagate with speed 
vp sa c/300, the effective coupling constant for elec- 
tromagnetic interactions turns out to be quite large, 
a = a /vp = 300/137 « 2 for suspended graphene 
[1—3]. In this case the non-perturbative effects could 
play an important role. Theoretical considerations sug- 
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gest that such strong interaction between charge carriers 
could result in the insulator-scmimetal phase transition 
in graphene [4-9]. 

However, due to the large value of the effective cou- 
pling constant there are no reliable analytical methods 
which allow to study this phase transition from the first 
principles, and one has to use numerical simulations. The 
effective field theory of graphene can be efficiently sim- 
ulated using lattice staggered fermions [10-17]. A single 
flavor of staggered fermions on 2 + 1-dimensional square 
lattice corresponds to two independent flavors of contin- 
uum Dirac fermions [18-20], which exactly reproduces 
the number of fermion flavors in the graphene effective 
field theory. 

In papers [10-17] the insulator-semimetal phase tran- 
sition was studied numerically by considering the 
fermionic "chiral" condensate ('iptp). Within the effec- 
tive field theory of Dirac quasiparticles, nonzero con- 
densate signals the opening of a gap in the quasiparticle 
spectrum, thus it plays the role of the order parameter 
for the semimctal-insulator phase transition. In [10—14] 
Coulomb interactions between fermions were modeled by 
a non-compact 3 + 1-dimensional Abelian lattice gauge 
field. It was found that the condensate is formed at a 
critical coupling constant of the non-compact gauge field 
j3 ~ 0.1. Motivated by the theoretical considerations of 
[9, 21], the authors of [15-17] have also studied a similar 
theory with a contact interaction instead of Coulomb po- 
tential. They have also found that fermionic condensate 
is formed at sufficiently strong coupling. 

However, the value of the conductivity has not yet been 
directly measured in numerical simulations. In this paper 
we report on direct numerical measurements of AC and 
DC conductivities of graphene within the effective field 
theory of Dirac quasiparticles. Our lattice regularization 
of this effective field theory is similar to the one used 
in [10-14] . The conductivity is obtained from the Green- 
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Kubo dispersion relations for the ground-state correlators 
of electromagnetic currents. These relations are inverted 
with the help of the Maximum Entropy Method [22, 23]. 

In agreement with the predictions of [4-17] we find that 
when a nonzero fermionic condensate is formed, the DC 
conductivity rapidly decreases. For the maximal value 
of the coupling constant used in our simulations, which 
corresponds to substrate dielectric permittivity e = 1.75, 
the DC conductivity is smaller than the DC conductivity 
in the weak-coupling limit by at least three orders of 
magnitude. 

The paper is organized as follows: in Section 2 we 
briefly review the effective field theory of graphene and 
discuss its lattice regularization, as well as suitable sim- 
ulation algorithms. In Section 3 we present and discuss 
our numerical results for the graphene conductivity and 
the fermionic condensate ( rpip ) . Section 4 contains some 
concluding remarks and the discussion of the obtained 
results. 



2. LATTICE REGULARIZATION OF THE 
EFFECTIVE FIELD THEORY OF GRAPHENE 



thus arrive at the following partition function: 

Z = J V$Vi)VA a exp {-1J d A x (cUo) 2 - 

-Jd 3 xi, f (r Q (d - igA ) - T t d)j V>/j , (3) 

where the effective coupling constant g 2 = e 2 /vp sa 
300/137 ~ 2. Finite temperature T can be introduced 
by imposing periodic boundary conditions in Euclidean 
time x° with the period |y. 

By virtue of the commutation relations [O a ,For.;] = 
with O a = 1, T 3 , T 5 , ir 3 r 5 the action of the effective field 
theory (3) has the global U(4) symmetry 



ipf -» exp [iO a ® r b a ab ) ijjf 



(4) 



where to — 1 and n - arc spin Pauli matrices which act 
on flavor index /. Finite temperature and chemical po- 
tential do not break this global symmetry on the level 
of the Lagrangian, however, it might be broken sponta- 
neously due to sufficiently strong Coulomb interactions 
[4-17]. 



2.1. Basic definitions 



2.2. Lattice action 



We start from the Euclidean path integral representa- 
tion of the partition function of the effective field theory 
of graphene [1-3] : 

Z = J V^fV^fVA^ exp(~\ J d A x {d^A u] f- 
d^x^f r (d - ieA ) </>/ - 



E 



d 3 x^ f T l (di-ievpAi) V>/ , (1) 



where A^ , fi = . . . 3 is the vector potential of the 3 + 1 
electromagnetic field, are Euclidean gamma-matrices 
and tpf, f — 1,2 are two flavours of Dirac fcrmions 
which correspond to the two spin components of the non- 
relativistic electrons in graphene. We have also taken 
into account that Dirac fermions propagate in 2 + 1- 
dimensional subspace at x 3 = with speed vf ~ 1/300. 

After rescaling of the coordinates and the vector po- 
tential 



x°/v F , A° Ai -> 



1 



A, 



(2) 



we conclude that the fluctuations of the spatial compo- 
nents Ai of the vector potential are suppressed by a factor 
1/vjp and we can set Ai = in practical calculations. We 



Following [10-17] we use staggered fermions [19, 20] 
in order to discretize the fermionic part of the action in 
(3). One flavor of staggered fermions in 2 + 1 dimensions 
corresponds to two flavors of continuum Dirac fermions 
[18-20], which makes them especially suitable for simu- 
lations of the graphene effective field theory. 

The action for staggered fermions coupled to Abelian 
lattice gauge field is 

x V=0,l,2 



E 

^=0,1,2 



^xOix 



(5) 



where the lattice coordinates x take integer values 



0...L b 



1 and x 3 is restricted to x 



0, f a 



is a single-component Grassman-valued field, a J;(1 = 
(— 1) :c ii+---+ x m-i j and 6 Xtfl are the link variables which 
are the lattice counterpart of the vector potential A^ (x) . 
For further convenience, we have also introduced the ma- 
trix elements D X:V of the staggered Dirac operator. The 
fields ^> x , ^ x satisfy periodic boundary conditions in spa- 
tial directions and anti-periodic boundary conditions in 
the Euclidean time direction. To account for the latter, 
we make a shift 6 x ,o ~^ &x,o + t in (5) at the time slice 
with x° = 0. 
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In order to recover the original spinor and flavor indices 
of the continuum Dirac fermions in (3), the lattice should 
be subdivided into the cubic blocks consisting of 2 x 2 x 
2 elementary lattice cells. Now the coordinates of all 
lattice sites can be enumerated as x^ = 2y u + rj^ , where 
= 0>1- We define the new fields on the lattice of y 
coordinates [11, 19, 20]: 



*,]/ = 



1 



4^2 



(6) 



where W VyV is the product of e along the the path 
which connects lattice sites with coordinates 2y and 
2y + rj, a = 1,2,3,4 is the Dirac spinor index and 
/ = 1, ...,4 is the flavor index. It can be shown that in 
terms of these new fields defined on the lattice with dou- 
ble lattice spacing the staggered fermion action (5) repro- 
duces the naive discretization of the continuum fermionic 
action in (3). However, there are additional terms which 
explicitly break the global U (4) symmetry of the contin- 
uum action in (3) down to its U (1) ® U (1) subgroup 
and which decouple only in the long-wavelength limit 
[11, 19, 20]. 

Since the fermion action (5) is restricted to the 2 + 1- 
dimensional subspacc with x — 0, not all components of 
[& y ]" are independent. Due to the absence of T 3 in (6) 
it satisfies the constraint 



r 3 r 5 %t 5 t 3 



(7) 



It is easy to check that in a representation of Euclidean 
gamma- matrices with I^Ts = diag (1, 1, — 1, — 1) this 
constraint implies the following block-diagonal form of 
the matrices 



A 
B 



(8) 



which is equivalent to two flavors of 4-component Dirac 
spinors. 

Now let us consider lattice discretization of the action 
of the electromagnetic field in (3). There exist two basic 
formulations of the U (1) lattice gauge theory: compact 
and non-compact. In order to exclude non-physical con- 
fining phase of the compact U (1) gauge theory [24] here 
we use the non-compact action for the gauge fields: 



S \&x 



x+i, 



( = 1 



(9) 



where summation over x now goes over the whole four- 
dimensional lattice. As discussed above, the fluctuations 
of the spatial components of the vector potential At (x) 
are suppressed in the effective field theory of graphene 
(3). Correspondingly, we also set to zero the spatial link 
variables 9 X 

In continuous space, the inverse lattice coupling con- 
stant /3 is related to the substrate dielectric permittivity 



= 1 = -^! 

g 2 Ane 2 



(10) 



where the factor ^p- takes into account the screening of 
the electrostatic interactions by the substrate. However, 
this relation can be modified due to finite lattice spacing 
effects such as the flavor symmetry breaking for staggered 
fermions [25, 26]. Generally, such effects tend to shift the 
phase transition towards the weak-coupling region [25, 
26]. We leave the study of such finite-spacing artifacts 
for future work. 

We note also that although the gauge field action (9) is 
non-compact, the fermionic action (5) is still "compact", 
that is, periodic in the variables 9^ (x). In general, it is 
impossible to couple the gauge field to lattice fermions in 
a non-compact way while preserving the gauge invariance 
of the theory. 

Since the fermionic action (5) is bilinear in the fermion 
fields, they can be integrated out in the partition function 
(3): 



exp(-S g [6 x ,o] -Sy >,..M',..'y,.,/ ) = 

Ve x , det(D[6 X!0 ])exp(-S g [e Xi0 \). (11) 



Thus we deal with the effective action 

S eff [9 X , ] = S g [9 x , a ] - lndet (D [9 X , ]) (12) 

which includes the determinant det (D [9 X: J) of the stag- 
gered Dirac operator D XtV [9 XjiJ \ introduced in (5). 

2.3. Simulation algorithm 

We use the standard Hybrid Monte-Carlo Method for 
generation of the configurations of the field 9 Xt o with the 
statistical weight exp(— 5 e // [^,o]) [11, 19, 20]. 

In order to calculate the determinant of the staggered 
Dirac operator in (11), we take into account that in the 
basis of even and odd lattice sites (which are defined 
as lattice sites with even or odd sum of all coordinates 
x° + x 1 + x 2 ) it takes the form [11, 19, 20] 



D[0 w ,o] = 



m 



D. 



e o 

m 



(13) 



with D\ a = -D 
to 



The determinant of D is thus equal 



det (D) = dot 



Dt -Do 



(14) 



which is a manifestly positive quantity. Note that the 
operator m 2 + D\ D eo acts only on the subspace of even 
lattice sites. We use the ^-algorithm in our simulations 
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[19, 20], in which the determinant (14) is represented in 
terms of a Gaussian integral over the pseudo-fermion field 

<j>x- 

det (m 2 + D\ Q D eo ) = J V4> x Vcf> x 
cxp f - ]T + D l D eo )~l <t> y J , (15) 

\ x,y / 

where the sum over x, y goes only over even lattice 
sites. The field <j> x is then stochastically sampled with 
the weight (15). To this end we generate the ran- 
dom field £ x according to the Gaussian distribution 

P [£,x] ~ exp ^— £x £x^J and then calculate cf> x = 

( m2 + D\ Deo) 1 £, y at the beginning of each Molec- 

v x ' v 

ular Dynamics trajectory [19, 20]. Nonzero mass term 
in (5) and (13) is necessary in order to ensure the in- 
vertibility of the staggered Dirac operator. Numerical 
results in the physical limit of zero mass were obtained 
by performing simulations at several nonzero values of m 
and by extrapolating the expectation values of physical 
observables to m — > 0. 

In order to speed up the simulations we also perform 
local heatbath updates of the gauge field outside of the 
graphene plane (at a; 3 ^ 0) between Hybrid Monte-Carlo 
updates. Both algorithms satisfy the detailed balance 
condition for the weight (11) [19, 20]. Successive ap- 
plication of these algorithms does not, in general, have 
this property. Nevertheless, by using the composition 
rule for transition probabilities it is easy to demonstrate 
that the path integral weight (11) is still the stationary 
probability distribution for such a combination of both 
algorithms. While local heatbath updates are computa- 
tionally very cheap, they significantly decrease the auto- 
correlation time of the algorithm. 

2.4. Physical observables on the lattice 

The main goal of this paper is to measure the elec- 
tric conductivity of graphene, that is, a linear response 
of the electric current density Ji (x) — tJj (x) 7, ip (x) 
to the applied homogeneous electric field Ej (t) (where 
t is the real Minkowski time). It is convenient to in- 
troduce the AC conductivity cr^ (w), so that Ji(w) — 
(Tij (w) Ej (w), where Ji(w) = j dt e~ vwt Ji (t) and 
Ej (uu) = J dt e~ twt Ej (t). Due to rotational symmetry 
of the effective field theory (3), Oij (w) should have the 
form &ij (w) = Sij a (w). Correspondingly, the DC con- 
ductivity is equal to the value of a (w) at w — > 0. 

By virtue of the Green-Kubo dispersion relations [22, 
23, 27], the Euclidean current-current correlators 

G{r) = \Y, I dxl dx2 ( J * (°) J * ( 16 ) 

i=l,2 



can be expressed in terms of a (w) as 

G(r) = J ^K(w,r)a(w), (17) 


where the thermal kernel K (w, r) is 

WCOSh (w (t — TPr)) , » 

K(w - r, = s ,lh( # ) (18 » 

and r = x° is the Euclidean time. We use here a non- 
standard definition of the kernel (18) from [23], which is 
more convenient for numerical analysis. 

Note that the current density in graphene is the charge 
which flows through the unit length in unit time and thus 
has the dimensionality of L~ 2 (where L stands for length) 
in units with H = c = 1. Correspondingly, the current 
density in lattice units is the charge which flows through 
a link of the dual lattice of length a in time cl/vf- Thus 
in order to express the current-current correlator (16) in 
physical units, one should multiply the result obtained 
on the lattice by a 2 v F /a 4 , where an additional factor of 
a comes from integration over x 1 , x 2 in (16). With the 
Euclidean time r in (16), (17) and (18) being expressed 
in units of lattice spacing in temporal direction cl/vf, 
integration over w in (17) also includes a factor v F /a 2 . 
We thus conclude that the AC conductivity a (w) is di- 
mensionless. Moreover, the DC conductivity a (0) is a 
universal quantity which does not depend on the lattice 
spacing or on the ratio of lattice spacings in temporal 
and spatial directions. For conversion to the SI system 
of units, it should be multiplied by e 2 / (2nh). 

In numerical simulations G (r) is measured for several 
(~ 10 1 ) discrete values of r. A commonly used method 
to invert the relation (17) and to extract the continuum 
function a (w) from the lattice discretization of G (r) is 
the Maximum Entropy Method [22, 23]. 

For staggered fermions the electric current Ji (y) can 
be expressed in terms of the fields ^ x as [20] : 

»? 

+*2 tf +,44°***a»-H/)> (19) 

where we have taken into account that the spatial link 
variables 9 x i are effectively equal to zero. Since the cur- 
rent (19) is defined on the lattice with double lattice spac- 
ing, we calculate the Euclidean current-current correlator 
(16) only on time slices with even r. 

In order to make sure that we reproduce the results 
of [10-14], we have also calculated the fermionic chiral 
condensate. In terms of staggered fermions it can be 
written as 

<W = 8Z^ (20) 
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After the fermions in the partition function are inte- 
grated out, the current-current correlator (16) and the 
chiral condensate (20) can be expressed in terms of ex- 
pectation values of certain combinations of the staggered 
fermion propagator D~^ y [9 XtfJ ] with respect to the weight 
(11). We give the explicit expressions for these combina- 
tions in Appendix A. 



3. NUMERICAL RESULTS 

Using the algorithm described in Subsection 2.3, we 
have generated 400 statistically independent gauge field 
configurations on the 20 4 lattice for each point in the 
space of lattice parameters P and to. For each value of 
P in the range (3 = 0.05 . . . 0.025 (which corresponds to 
substrate dielectric permittivities e = 1.75... 12. 75 ac- 
cording to (10)) the measurements were performed at 
three different values of mass m = 0.01, 0.02, 0.03. For 
the smallest mass to = 0.005, for which the simula- 
tions are most expensive, /? took values in the range 
P = 0.05... 0.15 (e = 1.75... 7.25). In order to esti- 
mate the finite-volume effects, we have also generated 
100 gauge field configurations on the 28 4 lattice with 
p = 0.05 ... 0.21 and m = 0.01. 



the data to the limit m —> 0. The result of such extrap- 
olation is also shown on Fig. 1. It suggests that there 
is a critical value of ft in the range 0.08 < p c < 0.09 
(which corresponds to 3.4 < e c < 4.0) such that the con- 
densate is zero for p > /3 C . A fit of the extrapolated 
data of the form (tpip) = b{p c — /3) 7 for p < p c yields 
Pc = 0.091 ± 0.002(e c = 4.0 ± 0.1) and 7 = 1.0 ± 0.16. 
These values are in agreement with the results of [10, 11], 
where the critical inverse coupling constant p c and the 
critical index 7 were estimated as 0.071 > p c < 0.091 
and 7 i= 1. 
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FIG. 1: Fermionic condensate (tptp) as a function of inverse 
lattice coupling constants ft at different values of mass m and 
extrapolation to the limit m — s> 0. Solid line is the fit of the 
extrapolated data with the function (ipip} ~ (j3 c — PY' with 
p c = 0.0908 ± 0.0018 and 7 = 1.0 ± 0.16. 

To check our simulation algorithms and to make sure 
that we reproduce the results of [10, 11], we first con- 
sider the fermionic chiral condensate (ipxp). It is plotted 
on Fig. 1 for different values of the mass to. The con- 
densate rapidly decreases as the inverse lattice coupling 
constant p (or, equivalently, the substrate dielectric per- 
mittivity e (10)) is increased. For each value of ft, we 
have fitted the mass dependence of the condensate with 
a quadratic polynomial and used this fit to extrapolate 



FIG. 2: Current - current correlators (16) for m — 0.01. 
Solid lines are the fits obtained using the Maximum Entropy 
Method. 

Our measurements of the conductivity of graphene 
start from the calculation of the current-current correla- 
tors (16), which are shown on Fig. 2 for to = 0.01 and for 
different values of inverse coupling constant p. We note 
that the correlators decay significantly faster for smaller 
values of p. 

The AC conductivity a (w) is extracted from the cor- 
relators using the Maximum Entropy Method with four 
basis eigenfunctions of the kernel (18) and a constant 
model function <tq (w) — 0.1 [22, 23]. The profiles of 
a (vS) are plotted on Fig. 3 and the corresponding fits of 
the correlators arc shown on Fig. 2 with solid lines. 

It is important to note at this point that on Fig. 3 
the angular frequency w is given in lattice units. For 
qualitative comparison with experimental data one can 
assume that the lattice spacing a for the spatial direc- 
tions of the cubic lattice used in our simulations is com- 
parable with the lattice spacing a — 0.246 nm of the 
hexagonal lattice in graphene. After the rescaling (2) 
the discretization step for the Euclidean time r should 
be of order At ~ a/vp. The temperature T in (17) and 
(18) is then equal to kT = H/ (Lq At) - 0.1 eV, which 
is much smaller than the characteristic binding energy 
in graphene ~ 1 cV. Thus our simulation results should 
correspond to sufficiently low physical temperatures as 
compared to characteristic excitation energies. 
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FIG. 3: AC conductivity a (w) in (17) for m = 0.01 at /3 close to /3 C (on the left) and at /3 > j3 c (on the right). 



For the inverse coupling constant j3 below approxi- 
mately 0.08 (e < 3.4), a (w) has one very broad peak 
around w ~ 1.2, and the DC conductivity a (0) has some 
small nonzero value. As /3 increases towards /3 C , the sec- 
ond peak emerges at w — and both peaks become nar- 
rower and higher (see Fig. 3, left plot). The emergence 
of the second peak results in the rapid growth of the DC 
conductivity. At p > fi c , the two peaks continue to grow, 
and their widths become comparable to the temperature 
T = Lq 1 in lattice units (see Fig. 3, right plot). 

In order to understand such peak structure in the 
weak-coupling limit, remember that for free Dirac 
fcrmions the AC conductivity a (w) has a delta-function 
singularity at w = [30, 31], which is a manifestation of 
the absence of scattering of charge carriers. When the 
interactions are turned on, this peak is smeared, which 
results in a large but finite value of the DC conductivity 
a (0). The second peak practically does not move as the 
mass to is changed. We conjecture that this second peak 
corresponds to a saddle point in the dispersion relation 
of staggered fermions which is situated in the middle of a 
straight line which connects the two distinct Fermi points 
[28] . The position of this peak should thus depend on the 
lattice regularization of the effective field theory (3) and 
should correspond to the optical frequency range for real 
graphenc. 

The DC conductivity a (0) is shown on Fig. 4 as a func- 
tion of inverse coupling constant /? at different values of 
the mass to. We normalize the conductivity to a single 
spin component and a single Dirac point, thus on Fig. 4 
we plot a (0) /4 rather than a (0). a (0) quickly decreases 
as both to and /3 become smaller. One can also note two 
distinct discontinuities in the dependence of a (0) on (3. 
For example, at m = 0.01 the first discontinuity is situ- 
ated between /3 = 0.07 (e = 3) and (3 = 0.09 (e = 4), and 
the second one - between /3 = 0.12 (e = 5.6) and f3 = 0.13 
(e = 6.2). The position of the first discontinuity depends 
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P 



FIG. 4: DC conductivity per spin per valley a (0) /4 in units 
of e 2 /h as a function of inverse lattice coupling constants at 
different values of mass m. The result of extrapolation of the 
data to the limit m — > is plotted with black dots and solid 
lines. 



only weakly on the mass m and roughly corresponds to 
the critical inverse coupling constant /3 C obtained from 
the analysis of the chiral condensate (see Fig. 1). The sec- 
ond discontinuity shifts to smaller (3 and becomes some- 
what weaker as to decreases. Linear extrapolation of the 
position of this discontinuity to the limit of zero mass 
(see Fig. 5) suggests that at to = both discontinuities 
coincide. We also note that the profile of the AC con- 
ductivity a (w) practically does not change across this 
second discontinuity. Thus there seems to be a single 
phase transition in the chiral limit to — ¥ 0, in agreement 
with the results of [10-12]. The corresponding critical 
value of the inverse coupling constant can be estimated 
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to lie in the range 0.07 < f3 c < 0.09 (3 < e c < 4). 

For each value of (3 between the two discontinuities, 
we also perform the quadratic fit of the mass dependence 
of the conductivity and use it to extrapolate the data to 
m — > 0. This extrapolation is shown on Fig. 4 with black 
dots and solid lines. 

On Fig. 6 we compare the DC conductivity at m = 0.01 
on 20 4 and 28 4 lattices. The dependence of a (0) on the 
inverse coupling constant f3 is qualitatively the same for 
both lattices, in particular, the positions of the discon- 
tinuities practically coincide. However, the actual val- 
ues of the DC conductivity differ beyond the error bars, 
especially in the strong-coupling phase. This suggests 
that finite-volume and finite-temperature effects could 
be quite large for our lattice parameters. Indeed, quite 
large finite-temperature effects have been reported in a 
recent Monte-Carlo study of the tight-binding model on 
the hexagonal lattice [28] , where the values of lattice pa- 
rameters were quite close to those used in this work. We 
leave the detailed study of finite-temperature and finite- 
volume effects as a direction for further investigations. 




0.000 0.005 0.010 0.015 0.020 0.025 0.030 



FIG. 5: Linear extrapolation of the position of the second 
discontinuity of the DC conductivity to the limit m — > 0. 



CONCLUSIONS 



■ 28* lattice 
• 20* lattice 



FIG. 6: DC conductivity per spin per valley a (0) /4 in units 
of e 2 /h as a function of inverse lattice coupling constant /3 for 
20 4 and 28 4 lattices with m = 0.01. 



scmimetal phase transition in graphene. Interestingly, 
our estimate of e c is close to the dielectric permittivity of 
silicon dioxide tsi0 2 = 3.9, which is often used as a sub- 
strate for graphene. According to the data presented on 
Fig. 4, for the largest value of the coupling constant used 
in our simulations (/? = l/<? 2 = 0.05) the DC conductiv- 
ity turns out to be smaller than the DC conductivity in 
the semimetal phase (at (3 > (3 C ) by a factor of order of 
10 3 . 

Finally, we note that our value of the DC conductiv- 
ity in the semimetal phase is significantly larger than 
the conductivity of non-interacting quasiparticles in the 
ideal monolayer graphene (a (0) ~ Ae 2 /h) obtained in 
[2, 29, 30] using the Landauer approach. However, as 
was stressed in [30], in Landauer approach a finite value 
of the conductivity of free Dirac fermions is determined 
solely by scattering on the boundaries of the sample. In 
the absence of boundaries (for instance, on the lattice 
with torus topology) the AC conductivity a (w) has a 
delta-function singularity at w = 0, and thus the DC 
conductivity a (0) is formally infinite. In an interacting 
theory, this singularity is smeared out, and the DC con- 
ductivity takes some finite (but large) value. 



In this paper we have numerically studied the AC 
and DC conductivities of graphene by using lattice 
Monte-Carlo simulations with 2 + 1-dimensional stag- 
gered fermions which interact with 3+1-dimcnsional non- 
compact Abelian lattice gauge field. We have found that 
in a phase with spontaneously broken chiral symmetry 
(which corresponds to sublatticc symmetry of the orig- 
inal hexagonal lattice) the DC conductivity rapidly de- 
creases as the substrate dielectric permittivity becomes 
smaller. The estimates of the corresponding critical val- 
ues 0.07 < (3 C < 0.09 (3 < e c < 4) obtained both from the 
measurements of the chiral condensate and the conduc- 
tivity agree with each other and with the results of [10- 
12, 14]. This supports the existence of a single insulator- 
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computing Center of the Moscow State University. where S (x, y) = D 



Appendix A: Calculation of fermionic observables 

Here we give explicit expressions for the vacuum expec- 
tation values of fermionic observables used in our simu- 
lations. 

Fermionic chiral condensate corresponds to the diago- 
nal elements of the staggered fermion propagator: 

o L/V\ Li\ L/o — 

x 

where (...) on the right-hand side denotes averaging over 
lattice gauge field Xffi with the weight (11). We calculate 
this trace using the stochastic estimator [19, 20]. 

Current-current correlator is a sum of connected and 
disconnected parts: 

G (y°) = C (y°) - D (y°) . (A2) 

The connected contributions can be expressed in terms 
of staggered fermion propagator as [20]: 

Vi,V2 n,ri' 

x (( S (2y + n, i + ?/) 5 (r/, 2j/ + t + 77) ) + 
+ (5 (2y + 77, r/) 5 (i + r/, 2y + » + v) > + 
+(5(2j/ + i + r7,l + r7 / )5(j7 / ,2y + r ? )) + 
+ ( S (2y + i + v, rf) s(i + rf, 1y + r?) )) , (A3) 



For the calculation of the connected part of the corre- 
lator (A2) we take into account that the solution of the 
linear equation \y = D yjZ ip z with Xy = $x,y yields the 
staggered fermion propagator D~y for all y. 

The disconnected part takes the following form: 



3/1, V2 n,ri' 

x (s (2y + i, 2y + r?) + S (ly + tj, 2y + i + r?))). (A4) 



In practice the disconnected part of the correlator (A2) 
is calculated using stochastic estimators [19, 20], simi- 
larly to the chiral condensate. In our simulations we 
have found that the disconnected part of the correlator 
is much smaller and much noisier than the connected one. 
Therefore we have neglected it in our measurements of 
the conductivity. 
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